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POLE-BASED APPROXIMATION OF THE FERMI-DIRAC 

FUNCTION 

LIN LIN, JIANFENG LU, LEXING YING, AND WEINAN E 

Abstract. Two approaches for the efficient rational approximation of the 
Fermi-Dirac function are discussed: one uses the contour integral representa- 
tion and conformal mapping and the other is based on a version of the multipole 
representation of the Fermi-Dirac function that uses only simple poles. Both 
representations have logarithmic computational complexity. They are of great 
interest for electronic structure calculations. 



1. Introduction 

Given an effective one-particle Hamiltonian H , the inverse temperature (3 = 
and the chemical potential the finite temperature single-particle density 
matrix of the system is given by the Fermi operator 

(1) p = 2(1 + cMP{H - Ai)))-^ = 1 - tanh(|(/f - ^)) , 

where tanh is the hyperbolic tangent function. 

In the last decade or so, the development of accurate and numerically efficient 
representations of the Fermi operator has attracted a great deal of attention in the 
quest for linear scaling electronic structure methods based on effective one-electron 
Hamiltonians. These approaches have a numerical cost that scales linearly with A'', 
the number of electrons, and thus hold the promise of making electronic structure 
analysis of large systems feasible. Achieving linear scaling for realistic systems is 
very challenging. Formulations based on the Fermi operator are appealing since 
this operator gives directly the single particle density matrix without the need for 
diagonalizing the Hamiltonian. 

From a computational viewpoint, one main issue is that the right hand side of 
(ID) is an operator-valued function. To evaluate this function, we have to replace 
or approximate it by something which can be computed directly without diago- 
nalization. Obvious candidates are polynomial or rational approximations. Such 
an approach was first introduced by Baroni and Giannozzi [T| and Goedecker and 
co-workers [316] (see also the review article [4]). Several improvements have been 
made since then, for example, in [21 [3l [H [TOl [TTJ [Mj . These are put broadly under 
the umbrella of "Fermi operator expansion" (abbreviated as FOE) . 

From the viewpoint of efficiency, a major concern is the cost for representing 
the Fermi operator as a function of j3^E (for finite temperature) or AE/Eg (for 
gapped systems) where /3 is the inverse temperature, AE is the spectral width of 
the discretized Hamiltonian matrix and Eg is the spectrum gap of the Hamiltonian 
around the chemical potential. Consider a finite temperature gapless system for 
example, the cost of the FOE proposed by Goedecker et al scales as j3AE. The fast 
polynomial summation technique introduced by Head-Gordon et al [lOl [11] reduces 
the cost to (/3AJ?)^/^. The cost of the hybrid algorithm proposed by Parrinello et 
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al in a recent preprint [2] scales as {(3AEy^^ . The cost was brought down to loga- 
rithmic scaling \n{PAE) in [12] using a multipole representation of the Matsubara 
expansion of the Fcrmi-Dirac function. 

The purpose of this article is to introduce two alternative rational expansions of 
the Fermi-Dirac function that use only simple poles and have computational cost 
that scales logarithmically. The first strategy is to use the contour integral and 
conformal mapping idea proposed recently in [8]. This will be presented in the 
next section. The other strategy is to borrow ideas from [15] and use a version 
of multipole expansion [12j that only involves simple poles. This will be discussed 
in Section [3l Numerical examples illustrating the efficiency and accuracy of the 
representations are discussed in Section ID 



2. Rational expansions based on contour integral 

Our first approach is an adaptation of the ideas proposed recently in [8] based on 
contour integral representation and conformal mapping. Let us first briefly recall 
the main idea of [8j. Consider a function / that is analytic in C\(— oo,0] and an 
operator A with spectrum in [m,M] C M"*", one wants to evaluate f{A) using a 
rational expansion of / by discretizing the contour integral 

(2) f{A) = :^ I f{z){z-A)-Uz. 

The innovative technique in [8] was to construct a conformal map that maps the 
stripe S = [—K,K] x [0, A''] to the upper half (denoted as of the domain 
= C\((— oo, 0] U [to, M]). This special map from t g 5 to z G fl^ is given by 



(3) z = VmAlf— u sn(t) sn(t|fc), k 



Here sn{t) is one of the Jacobi elliptic functions and the numbers K and K' are 
complete elliptic integrals whose values are given by the condition that the map is 
from 5' to 0+. 

Applying the trapezoidal rule with Q equally spaced points in [—K + iK'/2, K + 
iK'/2), 



(4) t,^-K+ — + 2 '-' ^' , 1<J<Q, 



iK' {j-\)K 
2 Q 

we get the quadrature rule (denote Zj — zitj)) 



(5) MA) ^ -J^^^ |, /Ma_^L!^iM5£M, 

TrQfc ^ (fc i-sn(<j))2 

Here en and dn are the other two Jacobi elliptic functions in standard notation and 
the factor cn(ij) dn(tj)(/s^^ — sn{tj))^'^\^mM /k comes from the Jacobian of the 
function z(t). 

It is proved in [Sj that the convergence is exponential in the number of quadrature 
points Q and the exponent deteriorates only logarithmically as M/m oo: 

(6) |1/(A) - /q(A)|| = o(e-'Q/0°g(W™)+3)). 
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To adapt the idea to our setting with the Fermi-Dirac function or the hyperbohc 
tangent function, we face with two differences: First, the tanh function has singu- 
larities on the imaginary axis. Second, the operator we are considering, P{H — fi), 
has spectrum on both the negative and positive axis. 

2.1. Gapped case. We first consider the case when the Hamiltonian H has a gap 
in its spectrum around the chemical potential such that dist(/^, <j{H)) = Eg > 0. 
Physically, this will be the case when the system is an insulator. 

Let us consider f{z) = tanh('|z^/^) acting on the operator A = {H — fj,)^. 
Now, f{z) has singularities only on (— oo, 0] and the spectrum of A is contained in 
[El Eli where 

Em— max \E — ul. 

Eea(H) 

We note that obviously Em < Ai?. Hence we are back in the same scenario as 
considered in [8] except that we need to take care of different branches of the 
square root function when we apply the quadrature rule. 

More specifically, we construct the contour and quadrature points Zj in the 
z-plane using parameters m = Eg and M = Efj. Denote g{£^) = tanh(/?^/2), 

+ 1/2 

~ iZj , and B = H fi. The quadrature rule is then given by 
(7) gqiB) = — Im ' ^ ' ' 



where the factors in the denominator come from the Jacobian of the map from 
z to ^. The number of poles to be inverted is iVpoie = 2Q. After applying we 
have a similar error estimate for g{B) 

(8) \\g{B) - gQ{B)\\ = o(e— 'Q/(2i°g(SM/iJ,)+3))_ 

In Fig. [U a typical configuration of the quadrature points is shown. The x-axis 
is taken to be E' — fi. We see that in this case the contour consists of two loops, 
one around the spectrum below the chemical potential and the other around the 
spectrum above the chemical potential. 

Note further that as the temperature goes to zero, the Fermi-Dirac function 
converges to the step function: 

(9) m -~ 

Therefore, the contribution of the quadrature points on the right half plane 
(Re ^'^ > 0) is negligible when /3 is large. In particular, for the case of zero 
temperature, one may choose only the quadrature points on the left half plane. 
The quadrature formula we obtain then becomes 

The number of poles to be inverted is then A^poic = Q- 
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Figure 1 . A typical configuration of the poles on a two- loop con- 
tour. Q = 30, Eg = 0.2, Em = 4 and ;3 = 1000. The red line 
indicates the spectrum. The inset shows the poles close to the ori- 
gin. The X-axis is E — n with E the eigenvalue of H. The poles 
with negative imaginary parts arc not explicitly calculated. 



We show in Fig. [5] a typical configuration of the set of quadrature points. Only 
one loop is required compared with Fig. [TJ 



2.2. Gapless case. The more challenging case is when the spectrum of H does 
not have a gap, i.e., Eg = 0. Physically, this corresponds to the case of metallic 
systems. In this case, the construction discussed in the last subsection does not 
work. 

To overcome this problem, we note that the hyperbolic tangent function tanh(^z) 
is analytic except at poles (21 — l)7r//3i, ^ S Z on the imaginary axis. Therefore, we 
could construct a contour around the whole spectrum of H which passes through 
the imaginary axis on the upper half plane between the origin and n/jSi and also on 
the lower half plane between the origin and —Tr/fii. Thus, we will have a dumbbell 
shaped contour as shown in Fig. [3l 

To be more specific, let us first construct the contour and quadrature points Zj 
in the z-plane as in the last subsection using parameters m = tt^//?^ and Al = 
EIj + 7rV/32. Denote = ±{zj - tt^//]^/^, g = tanh(/3^/2) and B = H - ji. 
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Figure 2. A typical configuration of the poles for zero temper- 
ature (/? = oo). Q = 30, Eg = 0.2 and Em = 4. The red line 
indicates the spectrum. The inset zooms into the poles that is 
close to the origin. The x-axis is E — fi with E the eigenvalue 
of H. The poles with negative imaginary parts are not explicitly 
calculated. 



The quadrature rule takes the following form 
(11) 9QiB) 



-2KWmM 



TiQk 



Im 



^ .g(g+)(g+-ij)-^cnfe)dnfe) 



\3 = ^ 



E 



C(fc-i-sn(t,))2 



When apply the quadrature formula, the number of poles to be inverted is A^poio = 
2Q. Fig.Oshows a typical configuration of quadrature points for Q ~ 30. The map 
^(z) = (z — 7r^//3^)^/^ maps the circle in the z-plane to a dumbbell-shaped contour 
(put two branches together). 

Actually, what is done could be understood as follows. Similar to [5], we have 
constructed a map from the rectangular domain [— 3A', K] x [0, K'] to the upper 
half of the domain 

U^{z\lmz> 0}\([-^m,^m] Uz[7r//?,oo)). 

The map is carried out in three steps, shown in Fig. (4] The first two steps use the 
original map constructed in [8], however with extended domain of definition. First, 
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Figure 3. A typical configuration of the poles on a dumbbell- 
shaped contour. Q = 30, Eg = 0, = 4 and f3 = 1000. The 
inset zooms into the part close to the origin. The red line indicates 
the spectrum. The black crosses indicate the positions of the poles 
of tanh function on the imaginary axis. The poles with negative 
imaginary parts are not explicitly calculated. 



the Jacobi elliptic function 



(12) u = sn(t) = sn(t|fc), k 



y/M/m + 1 



maps the rectangular domain to the complex plane, with the ends mapping to 
[1, k~'^] and the middle vertical line -K + i[Q, K'] to [-fc"\ -1]. Then, the Mobius 
transformation 

(13) z = ^/mM 



maps the complex plane to itself in such a way that [— fc ^,—1] and [l,fc ^] are 
mapped to [0, to] and [M, oo], respectively. Finally, the shifted square root function 



(14) e=(2-m) 



1/2 



maps the complex plane to the upper-half plane (we choose the branch of the 
square root such that the lower-half plane is mapped to the second quadrant and 
the upper-half plane is mapped to the first quadrant), in such a way that [0,m] is 
sent to z[0, ^/m] and [M, oo) is sent to (— oo, —^/M — m] U M — m, oo). The map 
can be extended to a map from [— 7if, K] x [0, K'] to the whole J7, in this case, the 
z-plane becomes a double-covered Riemann surface with branch point at to. 
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Figure 4. The map from the rectangular domain [—3K,K] x 
[0, K'] to the upper-half of the domain U. The map is constructed 
in three steps: t — > u — > z — > f . The boundaries are shown in 
various colors and line styles. 



Since the function g is analytic in the domain U, the composite function g{t) = 
g{^{z{u{t)))) is analytic in the stripe in the <-plane, and therefore, the trapezoidal 
rule converges exponentially fast. Using a similar analysis that leads to it can 
be shown that 

(15) 115(B) -5Q(B)|HO(e-c^'3/i°s(/3i^-)), 

where C is a constant. 

We remark that the construction proposed in this subsection also applies to 
the gapped case. In practice, if the temperature is high (so that /3 is small) or 
the gap around the chemical potential is small (in particular, for gapless system), 
the contour passing through the imaginary axis will be favorable; otherwise, the 
construction in the last subsection will be more efficient. 



3. Rational approximations based on the multipole expansion 

Another strategy for obtaining an efficient rational approximation for the Fcrmi- 
Dirac function for finite temperature is based on the multipole expansion, proposed 
recently in [12] . Let us first recall the construction of the multipole representation. 

Using the Matsubara representation (pole expansion) of the Fcrmi-Dirac func- 
tion, the density matrix is given by 

(16) p = 1 - 4Re 



^/3(if-/i)-(2Z-l)7rz- 



The summation in (|16p can be seen as a summation of residues contributed from 
the poles {(2Z — l)7ri}, with I a positive integer, on the imaginary axis. This suggests 
looking for a multipole expansion of the contributions from the poles, as was done 
in the fast multipole method (FMM) [7]. To do so, we use a dyadic grouping of 
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the poles, in which the 71-th group contains terms from / = 2" ^ to ^ = 2" — 1, 
for a total of 2"^^ terms. We decompose the summation in ([T6| accordingly. Let 
X = ii{H - [l). Then 



(17) 



E 



1 



00 2"-l 



= 1 ^ ' n=l;=2"-i ^ ' n=l 



The basic idea is to combine the simple poles into a set of multipoles at ^ = 
where Z„ is taken as the midpoint of the interval [2"~^, 2" — 1] 



(18) 



3 • T- 



- 1 



Then the Sn term in the above equation can be written as 

2"-l 



E 



(19) 



00 

i yf- 



2(1 — ln)TTi 



V-^^ 1 / 2{l-ln)m 

X - (2/„ - l)7ri Vx-(2/„-l)7ri 

1 / 2{l-lr,)TTi 

_^x- {21 - l)7ri \x - {2ln - l)7ri 



0^ 



E 



i=2'' 



Using the fact that x is real, the second term in (|19p can be bounded by 

2"-l 



2"-l 



E 



1 



2{l — ln)Tri 



^^^Jx-{2l-l)m x-{2L-l)m " J(2/ - 1)^| 



1 



- 27r 3^ ■ 



2L 



1 



Therefore, if we approximate the sum Sn by the first P terms, the error decays 
exponentially fast with P: 



(20) 



p-i 



. X V- 1 2{l-ln)T:i Y 

, JT-i ^ - (2^" - l)'^* ^0 - (2Z„ - 1)^1 



< 



1 1 
2^3^' 



The above analysis is of course standard from the view point of the fast multipole 
method [7]. The overall philosophy is also similar: given a preset error tolerance, 
one selects the value of P, the number of terms to retain in 5„, according to ([20|) . 

Moreover, the remainder of the sum in (I16p from I = 7\/poio + 1 to 00 has an 
explicit expression 



(21) Re V ^ 

^ ' ^ 2x-{2l-\)i'K 



2ti 



Im ?/' A/poic 



where "0 is the digamma function ^{z) ~ V [z) /T{z). 
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In summary, we arrive at the following multipole representation for the Fermi 
operator [T2| : 



(22) 

p=l-4Rc ^ ^ 



Ng 2^1 

^ 1^ f3{H - /i) - {21 



p-i 



n=l 1=2" 



2{l — /„)7ri 



-//3(ff-/.)-(2/„-l) 



„ — l)7ri/ 



-Im V ( Mpoic + 2 + - A^) ) + 0{Ng/3 



Here A'^g is the number of groups in the multipole representation. Mpoio = 2^^^ — 1 
is the number of poles that are effectively represented in the original Matsubara 
representation. In practice, Ng simple poles are first calculated, and then the 
Nq{P — 1) multipoles can be constructed through matrix-matrix multiplication. 

A disadvantage of ((22l) is that one needs to multiply simple poles together to get 
the multipoles before extracting the diagonal of Fermi operator. This prevents us 
from being able to directly apply the fast algorithms for extracting the diagonal of 
an inverse matrix, such as the one proposed in [T31. Therefore, it is useful to find 
an expansion similar to (|^^ that uses only simple poles. As we mentioned earlier, 
the key idea in deriving ([22)) is to combine the poles in each group together to 
form multipoles as the distance between them and the real axis is large. However, 
if instead we want an expansion that involves only simple poles, it is natural to 
revisit the variants of FMM that only use simple poles, for example, the version 
introduced in [T3]. The basic idea there is to use a set of equivalent charges on a 
circle surrounding the poles in each group to reproduce the effective potential away 
from these poles. 

Specifically, take the group of poles from / = 2"~^ to Z = 2" — 1 for example. 
Consider a circle Bn with center c„ = (3 • 2"~^ — 2)7ri and radius r„ = 2"~^7r. It is 
clear that the circle Bn encloses the poles considered. Take P equally spaced points 
{^n,k}k=i Oil the circle Bn- Next, one needs to place equivalent charges {pn,k}k=i 
at these points such that the potential produced by these equivalent charges match 
with the potential produced by the poles inside Bn away from the circle. This can 
be done in several ways, for example, by matching the multipole expansion, by 
discretizing the potential on Bn generated by the poles, and so on. Here we follow 
the approach used in [TS] . 

We simply take a bigger concentric circle Bn outside Bn with radius i?„ = 2"7r 
and match the potential generated on Bn by the poles and by the equivalent charges 
on Bn- For this purpose, we solve for pn.k the equations 

(23) E^^- J u - y^^- 



Regularization techniques such as Tikhonov rcgularization arc required here since 
this is a first-kind Fredholm equation. 

One can also prove that similar to the original version of the multipole rep- 
resentation, the error in the potential produced by the equivalent charges decay 
exponentially in P, the details can be found in [15]. Putting these all together, we 
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can write down the following expansion of the Fermi-Dirac function 

No P 

- hm ^ (^Mpoie + i + ^PiH - fi)^ + 0{Ng/3''). 

The number of poles that are effectively represented in the original Matsubara 
representation is still A/poic = — 1. A^poio = NqP simple poles are now to be 
calculated in practice. 

The tail part can be approximated using a Chcbyshcv polynomial expansion. 
Similar to the analysis in |12| , it can be shown that the complexity of the expansion 
is 0{\ogf3AE). As wc pointed out earlier, the advantage of ((24)) over ((22l) is that 
only simple poles are involved in the formula. This is useful when combined with 
fast algorithms for extracting the diagonal of an inverse matrix [13j . 

Note that in (|22p and ([M)) . for 2"~^ < P there would be no savings if we use 
P terms in the expansion. They are written in this form just for simplicity. In 
practice the first P simple poles will be calculated separately and the multipole 
expansion will be used starting from the [P + l)-th term and the starting level is 
n = log2 P + I. We show in Fig. [5] a typical configuration of the set of poles in the 
multipole representation type algorithm. 
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Re 

Figure 5. A typical configuration of the poles in the multipole 
representation type algorithm. Mpoio = 512 and P = 16 is used 
in this figure. The poles with negative imaginary parts are not 
explicitly shown. The inset shows the first few poles. The first 16 
poles are calculated separately and the starting level is n = 5. 
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Table 1 . A'^poic and error of electronic density per electron with 
respect to various (3AE. The energy gap Eg « 0.01. The contour 
integral representation for gapped system at finite temperature is 
used for the calculation. The performance of the algorithm depends 
weakly on (3AE. 



4. Numerical results 

We test the algorithms described above using a two dimensional nearest neighbor 
tight binding model for the Hamiltonian. The matrix components of the Hamilton- 
ian can be written as (in atomic units). 



(25) Hi>j>-ij 




i ± 1,/ = j or i' = = i ± 1. 



The on-site potential energy Vij is chosen to be a uniform random number between 
and lO"'^. The domain size is 32 x 32 with periodic boundary condition. The 
chemical potential will be specified later. The accuracy is measured by the error 
of the electronic density profile per electron 

Tr IP - PI 

(26) Ap„i = -r^ ^. 

JVElcctron 

4.1. Contour integral representation: gapped case. The error of the contour 
integral representation is determined by A^poio- At finite temperature iVpoic = 2(5, 
while at zero temperature A^poic = Q, with Q being the quadrature points on one 
loop of the contour. The performance of the algorithm is studied by the minimum 
number of iVpoio such that Ap^d (the error in the electronic density per electron) 
is smaller than 10~^. For a given temperature, the chemical potential /x is set to 
satisfy 

(27) TrP = ATgi^etron. 

In our setup the energy gap Eg « 0.01 Hartree = 0.27 eV and Em ~ 4 Hartree. 
Therefore, this system can be regarded as a crude model for semiconductor with a 
small energy gap. The number of A'^poie and the error Ap^ei are shown in Table [1] 
with respect to PAE ranging between 4, 000 and up to 270, 000. Because of the 
existence of the finite energy gap, the performance is essentially independent of 
PAE, as is clearly shown in Table [TJ 

When the temperature is low and therefore when (3 is large, as discussed before 
the finite temperature result is well approximated by the zero temperature Fermi 
operator, i.e., the matrix sign function. In such case the quadrature formula is 
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given by (jlOp . Only the contour that encircles the spectrum lower than chemical 
potential is calculated, and A^poio = Q- 

In order to study the dependence of Ap^ci on the number of poles A^poic, we tune 
artificially the chemical potential to reduce the energy gap to 10"^ Hartree. Fig. [6] 
shows the exponential decay of Apioi with respect to A^poio- For example, in order 
to reach the 10~® error criterion, A^poio ~ 50 is sufficient. The increase in A^poic is 
very small compared to the large decrease of energy gap and this is consistent the 
logarithmic dependence of A'^poio on Eg given by (H]). 




10 20 30 40 50 60 70 

-^Pole 

Figure 6. The lin-log plot of the error of electronic density 
per electron with respect to A^poic- The energy gap Eg « 10~^. 
The contour integral representation for gapped system at zero- 
temperature is used for calculation. 

4.2. Contour integral representation: gapless case. For gapless systems such 
as metallic systems, our quadrature formula in (jlip exploits the effective gap on 
the imaginary axis due to finite temperature. In the following results the chemical 
potential is set artificially so that Eg = 0. Em ~ 4 Hatree and the error criterion is 
still 10^^ as in the gapped case. Table [2] reports the number of poles A^poic and the 
error Aprei with respect to (3/S.E ranging from 4, 000 up to 4 million. These results 
are further summarized in Fig. [7] to show the logarithmic dependence of A'poic on 
pAE, as predicted in the analysis of (fT5)) . 
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Table 2. A'poic and error of electronic density per electron 
with respect to various (3AE. Eg = 0. The contour integral repre- 
sentation for gapless system is used for the calculation. 




Figure 7. Log-lin plot of iVpoio with respect to (3AE. The contour 
integral representation for gapless system is used for the calcula- 
tion. 



4.3. Multipole representation. The approach p4)) based on the multipole rep- 
resentation has three parts of error: the finite-term multipole expansion, the finite- 
term Chebyshev expansion for the tail part, and the truncated matrix-matrix mul- 
tiplication in the Chebyshev expansion. 
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22 


4.76 X 10" 


-7 


16,832 


2,048 


128 


22 


4.84 X 10" 


-7 


33,664 


4, 096 


144 


22 


4.88 X 10" 


'7 


67,328 


8,192 


160 


22 


4.90 X 10" 


-7 


134,656 


16,384 


176 


22 


4.90 X 10" 


-7 


269,312 


32, 768 


192 


22 


6.98 X 10- 


-7 


538,624 


65, 536 


208 


22 


3.20 X 10- 


-6 


1,077,248 


131,072 


224 


22 


7.60 X 10- 


-6 



Table 3. The number of poles calculated iVpoie, the order of 
Chcbyshev expansion for the tail part A^cheb, and the error 
of electronic density per electron with respect to various f3AE. 
The number of poles excluded in the tail part Mpoio is chosen to 
be proportional to (3AE. 



The error from the multipolc expansion is well controlled by P in (|24|) . When 
P ~ 16, 1/3^ ~ 0(10^*). The number of groups Nq is usually no more than 
20, and therefore the error introduced by multipole expansion is around 0(10^^), 
which is much less than the error criterion 10~^. 

The number of terms in the Chebyshev expansion for the tail part iVchcb is 
O ( ) , with il/poic being the number of poles that are excluded in the tail part in 

the pole expansion. The truncation radius for the tail part is C'(cxp(— C-^y-)) . In 
order to reach a fixed target accuracy, we set Mpoic to be proportional to (3AE. Due 
to the fact that Mpoic ~ 2^^ w 2^p°'"/-f', this requires A'^poic to grow logarithmically 
with respect to PAE. 

The target accuracy for the Chebyshev expansion is set to be 10^' and the trun- 
cation radius for the tail is set to be 4 for the metallic system under consideration. 
For f3AE = 4208, Mpoic is set to be 512 so that the error is smaller than 10~^. For 
other cases, Mpoio scales linearly with (3AE. The lin-log plot in Fig. [5] shows the 
logarithmic dependence of iVpoic with respect to (3AE. For more detailed results, 
Table [3] measures Mpoie, A^poie, -^chob, and Apici for (iAE ranging from 4000 up 
to 1 million. For all cases, iVcheb is kept as a small constant. Note that the trun- 
cation radius is always set to be a small number 4, and this indicates the tail part 
is extremely localized in the multipole representation due to the effectively raised 
temperature. 

Table [3] indicates that the error exhibits some slight growth. We believe that it 
comes from the growth of the number of groups in the multipole representation (j24[) 
and also the extra log log dependence on j3AE (see [12] for details). When compared 
with the results reported in Table [H we see that for the current application to 
electronic structure, the contour integral representation outperforms the multipole 
representation in terms of both the accuracy and the number of poles used. 

5. Conclusion 

We propose two approaches for the expansion of Fermi operator: a rational 
approximation based on the contour integral idea introduced in [5] and a variant of 
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Figure 8. log-lin plot of TVpoie with respect to l3AE. The multi- 
pole representation is used for the calculation. 

the multipole representation in [12j using only simple poles. Both approximations 
result in logarithmic scaling complexity with respect to /?Ae with small prefactor. 
Fast algorithms for electronic structure calculations can be obtained by combining 
these approaches with the algorithm introduced in [13] for extracting the diagonal 
of the inverse of a matrix. 
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